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Abstract 

We describe results of the cluster algorithm Special Purpose Processor simulations 
of the 2D Ising model with impurity bonds. Use of large lattices, with the number 
of spins up to 10^, permitted to define critical region of temperatures, where both 
finite size corrections and corrections to scaling are small. High accuracy data 
unambiguously show increase of magnetization and magnetic susceptibility effective 
exponents (3 and 7, caused by impurities. The M and x singularities became more 
sharp, while the specific heat singularity is smoothed. The specific heat is found 
to be in a good agreement with Dotsenko-Dotsenko theoretical predictions in the 
whole critical range of temperatures. 



1 Introduction 



The problem of inhomogenities influence on phase transitions has a long history. 
The first model with some kind of impurities to be studied, is, due to its simplicity, 
the Ising model. Effect of randomness on the critical behaviour of different Ising 
models has been investigated by Harris |^. He found that if the specific heat of the 
pure system diverges as some power of (T^ — T)~^, then the critical behaviour will 
be changed by impurities. Such a result does not give any information for the 2D 
Ising model, which has a logarithmic divergence of the specific heat. 

Theoretical treatment of 2D Ising model with ferromagnetic impurity bonds was 
pioneered by Vl.Dotsenko and Vik.Dotsenko (DD) [2-6]. They predicted new critical 
behaviour of the specific heat 0, ^, spin-spin correlation function, magnetization 
and magnetic susceptibility 0, ^. Later the same model has been considered in a 
number of theoretical works [7-15]. Some authors claimed the specific heat to remain 
finite for all temperatures [0, |13[. The others [7-12] confirmed DD result for the 
specific heat, but stated the critical behaviour of spin-spin correlation function and 
magnetization to be the same, or almost the same, as in the pure case, only slightly 
changed by logarithmic corrections. 

Experiments on quasi-2D compounds with almost Ising spins [16-20] did not 
show any deviations from Onsager results, probably because of unsufficient accu- 
racy, caused by large scale inhomogeneities, which smooth out the phase transition. 
Moreover, in such experiments it is not possible to exclude 3D effects. So, to un- 
derstand, if any of the theories catch essential physics, it is necessary to perform 
computer simulations. 



Early Monte Carlo simulations [21, ^ demonstrated the same critical behaviour 
as in the pure case. 

The large-scale simulations were started by Andreichenko et al [23-25]. The 
specific heat at the critical point of the model with impurities C{Tc) was studied as 
a function of the system size L [^ ^ . For large impurity strength the specific heat 
C{Tc) was found to be proportional to loglog(L), which seems to be in agreement 
with DD prediction. 

The behaviour of magnetization M(Tc) and magnetic susceptibility x{Tc) as 



a function of L turned out ^ to be the same as in the pure case, which 
contradicted to DD theory. 

Usually critical behaviour is studied as a function of relative temperature distance 
r = {Tc—T)/Tc from Tc. Some deviations from the pure behaviour of M(r) and x{^) 
were found by Wang et al [53], and in Special Purpose Processor (SPP) simulations 



26| , |27[] . But the accuracy of [25-27] was not very high, and the critical region was 



not defined. 

The new SPP |2^, |2^ , using cluster Wolff algorithm, permits to get very accurate 
description of the 2D random Ising model critical behaviour. It gives us not only 
the thermodynamic quantities, but also the spin-spin correlation function [^, ^ 

CF{r) =< S{0)S{r) > , 

which is directly studied by theories [4-12]. 
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2 The model 



We study the same system which has been considered earher [25-27]. Ising spins are 
located at the nodes of two-dimensional square lattice. To avoid appearance of the 
border-induced terms, we use periodic boundary conditions in both directions. 

Disorder is introduced by random distribution of impurity coupling constants 
over lattice bonds. Pure case exchange interaction constant is denoted as J, and 
impurity coupling constant is J' . The probability to find on some bond J' is p, and 
the probability to find J is (1 — p). In this paper we do not consider spin glasses, 
so both J and J' are ferromagnetic. 

What can be said about the phase diagram of such a model ? At p = there 
is the Onsager phase transition at (1/Tc) ~ 0.44068, if the pure coupling constant 
J is chosen to be equal to 1. The shift of the due to impurities for small p was 
calculated exactly in [0, ^ for J' = 0, and in ^ for the general case. 

Certainly, impurities destroy thermal fluctuations and decrease T^. If there are 
enough of strong impurities, then the phase transition disappears. Indeed, if J' = 0, 
and p > 0.5, then there is no percolation of bound spins in the system |3^. The 



lattice is broken into finite islands of interacting spins, and in a finite system there 
can be no phase transition. 

Fortunately, for p = 0.5 and J' > there exist exact duality relation [^, 34 1, 
which gives the value of as a function of J and J': 

tanh( J/T,) = exp(-2J7Tj . (1) 



In [^, ^ it was checked, that indeed, the position of the specific heat maximum 
tends to this Tc in the limit L — > oo. This confirms that there is only one phase 
transition point in the system. 

In their theory of the impure Ising model DD introduced a small parameter 

g (xp{J - J'f 

and the impurity induced length /, 

log(Zi) oc - . 
9 

The influence of impurities should be important only on the distances larger than 
li. If the disorder is small, g ^ 1, then for the finite lattice with the linear size L 
it may happen that li ^ L. In this case it is impossible to notice the influence of 
impurities on the critical behaviour. 

For this reason to study the deviations caused by inhomogeneities, we must 
choose g as large as possible to decrease /j. So, in real MC simulations g cannot 
be very small, and exact theoretical formula for g by DD 0, is no more valid. 
Nevertheless, as we show later, the DD formula for the specific heat 

C(T)(x-log(l + <7log(i)), (2) 
9 ^ 

where g is regarded just as a parameter, reasonably describes the simulation data. 
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3 Simulation procedure 



We used for the simulations the first cluster algorithm SPP ^|. It implements 
in hardware very efficient cluster Wolff algorithm [35-37], which is the improved 
version of the Swendsen - Wang algorithm P3 . The Wolff algorithm does not suffer 



of the critical slowing down, and at the critical point it should be about times 
faster than the conventional one spin-flip algorithm. 

In the case of our SPP L can be as large as 1024. So, the improvement in speed 
is about one million times. 



Detailed description of the SPP structure and functioning can be found in [29 
The class of problems, which can be solved by the SPP, is described in [^8|, . 



Precise meaning of the critical slowing down absence was found in the simulations 
of the pure case [^ : the relaxation time, measured in real simulation time, is the 



same at the critical point and far from it. 

The relaxation time is defined in the following way: we start simulation with all 
spins pointing in one direction. After some time all the thermodynamic values, such 
as magnetization or the lattice energy, come into the zone of thermal fluctuations 
near thermal equilibrium. We call this time the relaxation time. 

For the pure case it is necessary to flip about 20 Wolff clusters to get to the 



fluctuation region near Tc, and about 5 Wolff clusters far from Tc for L = 1024 ^9 
On the other hand, mean number of spins in the cluster far from the critical point 
is almost equal to the total number of spins L^, and near the number of spins in 
the cluster is about 4 times lower. So, the relaxation time, measured as the really 
spend computer time, is the same. 

The situation in the disordered system can be seen on Fig.l and 2. All results 
are given for the case J = 1, J' = 0.25, p = 0.5. According to (0) this corresponds 
to (1/T,) = 0.80705186. 

Fig.l shows the relaxation of magnetization and the correlation function CF{L/2) 
near Tc, at (1/Tc) = 0.808. We see that again, like in the pure case, 20 clusters should 
be flipped to enter the thermal fluctuation zone near Tc. 

As can be seen from Fig.l, there is some correlation between M and CF{L/2). 

Indeed, there are two ways to deflne M. First, we can count the difference in 
the number of up and down A^j^ spins in the system, then the magnetization is 
given by 

M, = iN^--N^)/iN^ + N^). 

The second deflnition of M is usually used to flnd it theoretically for the infinite 
system: 

M= (< 5(0)5(oo) >)i/2 _ 
In the finite system we can alternatively define M as 

Ms = (CF(L/2))i/2 _ 

Our simulations show that Mi and M2, averaged over impurity distribution, nor- 
mally lead to the same mean values of magnetization. Noticeable difference between 
them appears only very close to Tc, when the finite lattice size effects come into play. 
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Nevertheless, the correlation between M and CF{L/2) persists up to T^, as can be 
seen from Fig.l. 

The relaxation time is about 20 clusters not only for M, but also for other 
thermodynamic quantities. That can be seen in Fig. 2, which describes the relaxation 
of the neighbour spins correlation function CF(1). In the pure case CF(1) coincides 
with the energy per one bond. Now it is not so, because the energy is given by the 
mean value of 

< JoiS{0)S{l) > 

and there are some correlations between the value of coupling constant Jqi on the 
bond, connecting two neighbour nodes and 1, and the sign of the two spin product. 
Nevertheless, the relaxation curve for the energy, also shown in Fig. 2, behaves very 
much like relaxation curve for CF{1). 

Fluctuations of the magnetization, shown on Fig.l, are very large, as should be 
the fluctuations of the order parameter near the critical point. On the other hand, 
the energy, which is not the order parameter, does not fluctuate so strongly. 

Because the relaxation time near in the impure case is the same as in the pure 
case, we again come to the conclusion that the critical slowing down is absent. 

Nevertheless, to get all the thermodynamic data, described below, we permitted 
spins of the each sample to relax during first 2000 cluster flips. Only after that 
measurements for each sample were started. 

Each sample has its own distribution of impurity bonds. Coupling strengths 
J = 1 and J' = 0.25 were ascribed to each bond with the probability one half to 
insure selfduality condition. It is the difference between different samples which 
determines the value of standard errors for all the thermodynamic data. 

We use 3 sets of data. Low accuracy data, obtained with 10 samples, describe 
large r region. Better data, obtained with 100 samples, give general picture for all 
temperatures. Finally, very accurate 1000 samples data were obtained in the critical 
region, where asymptotics can be defined. 

The importance of the proper determination of the critical region can be seen 
from Fig. 3 and 4, on which we show reduced magnetization and magnetic suscep- 
tibility for the pure case, L = 1024 |^|. In the pure case the critical behaviour is 



described by power laws 

Mo = 1.22241 (r)^/^ (3) 

Xo = (0.025537 - 0.001989 r') {r'y^^^/T (4) 

where r' = {T,-T)/T 

Fig. 3 and 4 show ratios of the pure case ^9) M and x to this asymptotic laws. 



It is clear, that the critical region, in which the asymptotics are valid, is not very 
wide: 

0.001 <T< 0.02 (5) 
Low r restriction comes from finite lattice effects, which become important for 

r < (1/L) . (6) 
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On the other hand, we know the exact solution for the magnetization of the infinite 
system 

°° ^ sinh^(2/T)^ ■ 

This solution is true at large r, and shows that there should be analytic corrections 
to simple scaling law (H). These corrections to scaling lead to the large r limit in 



4 Results 

In the impure case asymptotics are no longer valid. Nevertheless, deviations 
from them are not very large. To see these deviations better, it is convenient again 
to divide M and x by Mq and xo- Corresponding ratios are shown as a function of 
T on Fig. 5 and 6. 

From this figures it is obvious, that the critical behaviour is changed, and can 
be described by larger effective exponents than in the pure case. It is also obvious, 
that the critical region for the impure case is somewhere in the limits 

0.003 <T < 0.03 . (7) 

To check once more, that small r behaviour is determined by the finite lattice size, 
we show the change of this behaviour with L on Fig. 7 and 8. From Fig. 8 we see that 
the maximum of x is shifted from r = 0.003 for L = 1024 to r = 0.006 for L = 512, 
as can be expected. 

It is natural to investigate critical region (|^) more carefully. The M and x 
results, obtained in this region using 1000 samples for L=1024, are shown in Fig. 9 
and 10. They can be described by changed effective exponents in power laws (0,^. 
This does not contradict to the description in terms of logarithmic corrections to 



that power laws p5|, but requires less number of fitting parameters. The effective 
exponents should not be regarded as real exponents, because that will violate the 
standard scaling relations. 

Fig. 11 shows the same data, as in Fig.9, additionally divided by r^, for e = 0.009, 
0.0075, 0.006. We see, that the effective exponent of M in the critical region is 
increased by 0.0075 from the pure Ising value 0.125. 

Fig. 12 shows the analogues data for the magnetic susceptibility x- In this case 
e = —0.11, —.135, —.17. This implies that the effective critical exponent for x is 
1.75 + 0.135. 

The specific heat C(r) is more difficult to study than the magnetization for two 
reasons. 

1) The specific heat is obtained as fluctuations of the energy according to the 
formula 

C =< {E~ <E>f> /T^ . 

As a result, fluctuations of C(r) for a given r do not decrease with the increasing 
lattice size L, as do the fluctuations of M. In reality, standard deviations of C(r) 
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depend only on r and the number of flipped clusters, which was used to measure 
C(r), but not on L. 

For this reason it has a sense to get large r values of C(r), for which the influence 
of finite lattice size is not important, on smaller lattices, which requires less computer 
time. 

2) The specific heat r dependence is complicated. This creates difficulties in 
defining the critical region. 

Despite these two problems, the DD formula (||) can help us to interpret the 
simulation data. 

Fig. 13 describes general behaviour of the specific heat in both pure and impure 
cases as a function of r. We see, that impurities reduce C(r) and cause deviations 
from the simple log(r) asymptotic behaviour. 

It is convenient to study the difference between C(r) and supposed asymptotics 
of it. 

Pure case asymptotic of C(r) per one node for r ^ is 

CP«-(r) = 0.4945 log(-) (8) 

T 

plus some constant. Analytic "corrections to scaling" can be made visible, if we draw 
the difference Upurei^) between the exact solution for the infinite lattice C^xacti''') 
(^. This difference 

Vpureir) = C'Z^tir) - 0.4945 log(i) + 0.6 (9) 

r 

is shown in Fig. 14. Small deviations from horizontal line at large r demonstrate 
"corrections to scaling" for the simple logarithmic behaviour of C^xacti'^)- 
The impure data in the same Fig. 14 are described by the curve 

y(r)=C(r)-0.211og(i) (10) 

T 

The coefficient 0.21 before log(r) in this formula was chosen to make ?/(r) curve 
horizontal near r = 0.01. We see, that deviations from the pure log behaviour are 
larger in the impure case. 

For this reason we may try to approximate C(r) by (@). It is natural to try to 
choose proportionality coefficient in such a way that for g ^ this formula turns to 
the pure Ising formula. Then C(r) must be chosen as 

4945 1 

C(r) = log(l + ^log(-)) + const , (11) 

9 T 

where g and const are parameters to be found from comparison with the simulation 
data. 

Fig. 15 shows 

4945 1 
<r) = C{t) - log(l + 0.295 log(-)) + 0.6 . (12) 
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Here g = 0.295 was chosen to make the curve as horizontal as possible in the critical 
region r > 0.003. Comparison with the pure case curve on the same Fig. 15 shows 
that "corrections to scaling" for the DD formula ([TI|) are approximately of the same 
value as in the pure case. 

The choice of g = 0.295 is justified by Fig. 16, which shows 1000 samples data for 
L = 1024 in the critical region. This Fig. 16 gives curves z{t) for 3 different values 
of ^ : 0.31, 0.295, 0.28. 

We can conclude that DD formula ([TTl) is surprisingly good in description of the 
impure case specific heat. The value of (7 is in a reasonable agreement with that 



found from logarithmic correction fit of magnetization data in [25 



5 Conclusions 

More pronounced singularities of magnetic properties M(r) and x(r) in the impure 
case are in qualitative agreement with the statement that x and M should remain 
the same functions of the correlation length ^, as in the pure case 



X oc 



f , (13) 



(14) 

but ^ grows faster than (l/r) when r ^ 0. 

So, the impurities increase the correlation length ^(r) for a given r. At a first 
glance, that seems to contradict to common sense. But we should take into account, 
that the impurities decrease the critical temperature T^, so the increased ^(r) is 
obtained at much lower absolute temperature T. 

It would be natural to say that our data give increased values of magnetization 
critical exponent P and magnetic susceptibility critical exponent 7. But this would 
violate the standard scaling formula a + 2/5 + 7 = 2. For this reason we regard 
increased f3 and 7 as effective exponents. This effective exponents can be treated as 
describing logarithmic corrections to the correlation length ^ ^ (1 — (7log(r))^/^/r 
P], combined with large analytic in r corrections to the scaling laws (p!3|,p^) . 
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Figure captions 



Fig.l: Time relaxation of ttie magnetization M (solid triangles) and of the corre- 
lation function CF{L/2) (empty triangles) as a function of a number Nc of flipped 
Wolff clusters. Broken lines show mean values, obtained for a large number of flipped 
clusters for the linear lattice size L = 1024. Temperature T is very close to Tc. 

Fig. 2: Time relaxation of the neighbour spin-spin correlation function CF(1) and 
of the energy per one bond (insert). Broken lines have the same meaning as in Fig.l. 



Fig. 3: Ratio of the magnetization M, obtained by the cluster SPP p9| for the pure 
Ising model, L = 1024, to the asymptotic law Mo(r). 



Fig. 4: Ratio of the magnetic susceptibility x, obtained by the cluster SPP for 
the pure Ising model, L = 1024, to the asymptotic law Xo{''~')- 

Fig. 5: Ratio of the random bond Ising Model magnetization to the pure case asymp- 
totic law Mq{t). Down triangles show 10 samples data, up triangles - 100 samples 
data. 

Fig. 6: Ratio of the random bond Ising Model magnetic susceptibility to the pure 
case asymptotic law Xo(t')- Down triangles show 10 samples data, up triangles - 
100 samples data. 

Fig. 7: Impure case M(t)/Mq{t) for two different lattice sizes. Empty circles show 
L = 512 data, solid diamonds - L = 1024 data. 

Fig. 8: Impure case x{'t') / Xoi'^') for L = 512 (empty circles) and L = 1024 (solid 
diamonds) . 

Fig. 9: 1000 samples critical region data for the impure case M{t)/Mq{t), L = 1024. 

Fig. 10: 1000 samples critical region data for the impure case x(t')/Xo(t')) L = 
1024. 

Fig. 11: Impure case M/(Mqt'^) for e = 0.006 (down triangles), 0.0075 (diamonds), 
0.009 (up triangles). 

Fig. 12: Impure case x/iXoi'^'Y) for e = —0.11 (down triangles), —0.135 (diamonds), 
—0.17 (up triangles). 
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Fig. 13: Specific lieat per one node C(r) for tlie L = 1024 impure system (up trian- 
gles), L = 512 impure system (empty boxes), and infinite pure Ising system (broken 
line) . 

Fig.l4: Difference y(r) between C(r) and tlie best logaritfimic approximation for 
it. Pure case ypure{T) is shown by the broken fine. Up triangles describe L — 1024 
impure data, empty boxes - L = 512 impure data. Number of samples is equal to 
100 in both cases. Only the dependence or independence of y on r is important, so 
all the differences in Fig.l3-Fig.l6 are displaced by arbitrary constant. 

Fig. 15: Difference z{t) between C(t) and the DD approximation for g — 0.295. Up 
triangles describe L = 1024 data, empty boxes - L = 512 data. Broken line shows 
infinite pure system C(r) deviations from the simple log(r) low. This deviations are 
caused by analytic corrections to scaling. Decrease of z{t) for r < 210~^ is caused 
by finite lattice effects. 

Fig.l6: 1000 samples, L — 1024, critical region difference z{t) between C(r) and 
DD approximations for g — 0.28 (down triangles), 0.295 (diamonds), 0.31 (up tri- 
angles) 
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